% Obtain bootstrapped standard errors
% According to settings in PAR structure, including the sample (main or
% high), the number of bootstrap iterations, and details about the
% estimation method
% 
% Called from the Main script.
% Input:  data, Results (estimation results / point estimates), and PAR
% Output: BS structure with std. errors in the field BS.BS_SE 
% BS saved to Output/Bootstrap_results_[main/high].m unless "outputfolder"
% is updated. 
% 
% Can be parallellized across compute nodes
function[BSResults,BSidx,noblow,BS,BSResultsorig] = bootstrap_func(data,PAR,Results,outputfolder)

    PAR.optsv0star.Display = 'off';
    PAR.bootstrap =1;
    PAR.optsv0star.MaxFunctionEvaluations = 50000;
    
    % Sample from the main or high-end sample only, as relevant
    if PAR.highvalue == 0
        weights = data.Sample_main;
        Nobs = sum(data.Sample_main);%to avoid oversampling
        sample="main";%for output file
    else
        weights = data.Sample_highend;
        Nobs = sum(data.Sample_highend);%to avoid oversampling
        sample="high";%for output file
    end

    [~,BSidx]=bootstrp(PAR.nrbootstraps,[],data,'Weights',weights);
    BSResultsorig = Results;
    parfor b = 1:PAR.nrbootstraps
        BSResultsorig(b)=estimation_func(data(BSidx(1:Nobs,b),:),PAR);
    end

    % Get standard errors and CI interval - report in table "BS"
    % Some models "blow up" - evidenced by huge or 0 std devs
    params       = fieldnames(BSResultsorig);
    npar         = numel(params)-1;%
    params=params(1:npar);
    BS_95CIlow   = NaN(npar,1);
    BS_95CIhigh  = NaN(npar,1);
    BS_SD        = NaN(npar,1);
    BS_SE        = NaN(npar,1);
    Res          = NaN(npar,1);
    
    % Exclude unstable samples (especially in limited high value sample)
    noblow       = [BSResultsorig.mus]<999;  
    BSResults    = BSResultsorig(noblow);
    N=sum(noblow);

    % Obtain statistics
    for i = 1:npar
        try
          par = [BSResults.(params{i})];
          BS_SD(i)       = std(par);
          BS_SE(i)       = BS_SD(i)/sqrt(N);
          BS_95CIlow(i)  = quantile(par,0.025);
          BS_95CIhigh(i) = quantile(par,0.975);
          Res(i)         = Results.(params{i});
        catch
          BS_SD(i)       = NaN;
          BS_SE(i)       = NaN;
          BS_95CIlow(i)  = NaN;
          BS_95CIhigh(i) = NaN;
          Res(i)         = Results.(params{i});
        end
    end
    BS = table(Res,BS_SD,BS_SE,BS_95CIlow,BS_95CIhigh,'RowNames',params);   
    
    % Save to output folder
    save(strcat(outputfolder,'Bootstrap_results_',sample),"BS")
            
   end
    